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Recent experimental evidence suggests that acute myeloid leukaemias may 
originate from multiple clones of malignant cells. Nevertheless, it is not known 
how the observed clones may differ with respect to cell properties, such as pro- 
liferation and self-renewal. There are scarcely any data on how these cell 
properties change due to chemotherapy and relapse. We propose a new math- 
ematical model to investigate the impact of cell properties on the multi-clonal 
composition of leukaemias. Model results imply that enhanced self-renewal 
may be a key mechanism in the clonal selection process. Simulations suggest 
that fast proliferating and highly self-renewing cells dominate at primary diagno- 
sis, while relapse following therapy-induced remission is triggered mostly 
by highly self-renewing but slowly proliferating cells. Comparison of simulation 
results to patient data demonstrates that the proposed model is consistent with 
clinically observed dynamics based on a clonal selection process. 
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1. Introduction 

Leukaemia is a clonal disease of the haematopoietic system leading to extensive 
expansion of malignant cells that are non-functional and cause impairment of 
blood cell formation. Recent experimental evidence indicates that the malignant 
cell population might be composed of multiple clones [1], maintained by cells 
with stem-like properties [2,3]. A clone consists of genetically identical stem and 
non-stem cells. Relapse of the disease after therapy is a common problem of 
leukaemias [1]. 

To understand better the origins of acute leukaemia relapses, a genetic inter- 
dependence between clones at diagnosis and relapse has been investigated 
using gene sequencing and other techniques. In most cases of acute lympho- 
blastic leukaemia (ALL), the clones dominating relapse were already present 
at diagnosis but were undetectable by routine methods [4-6]. Owing to quies- 
cence, very slow cycling or other intrinsic mechanisms [5,6], these clones 
survive chemotherapy and eventually expand [5,6]. This implies that the 
main mechanism of relapse in ALL is based on a selection of existing clones 
and not an acquisition of therapy-specific mutations [5]. Similar mechanisms 
have been described for acute myeloid leukaemia (AML), where clones at 
relapse are genetically closely related to clones at primary diagnosis [1,7] and 
did not have to acquire additional mutations during the course of disease [8,9]. 

Based on these findings, the evolution of malignant neoplasms can be inter- 
preted as a selection process [10-12] of cells with properties that enable them 
to survive treatment and to expand efficiently. Cells with different mutations 
may have different growth properties [1]. Chemotherapy significantly alters 
growth conditions of cells and therefore, it may have a strong impact on the 
selection process. If cells dominating at diagnosis are sensitive to therapy, 
minor clones with intrinsic resistance [5,6,13] may expand more efficiently 
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Figure 1. Schematic of the models. In Model 1, [a) self-renewal of leukaemic and haematopoietic cells depends on the total number of post-mitotic cells via 
negative feedback. In Model 2, {b) self-renewal of haematopoietic cells depends on mature blood cell counts. Leukaemic cells are independent of haematopoietic 
feedback signals. Increasing cell numbers in bone marrow space lead to increasing death rates of all bone marrow cell types, namely mitotic haematopoietic and all 
leukaemic cells. (Online version in colour.) 



once the competing clones are eliminated by the treatment. 
The latter could explain manifestation of different cell 
clones at diagnosis and at relapse without a need for 
additional mutations in between. 

The mechanism of the underlying selection process and its 
impacts on the disease dynamics and on the response of cancer 
cells to chemotherapy are not understood. Gene sequencing 
studies allow the genetic relationship between different clones 
to be deciphered, nevertheless the impact of many detected 
mutations on cell behaviour remains unclear [1], and often 
passenger mutations cannot be distinguished from relevant 
genetic changes [4]. Many authors, e.g. [14,15], have provided 
evidence for the heterogeneity of leukaemia stem cells (LSCs) 
attempting to identify LSC characteristics, for review see [16]. 
This heterogeneity is further supported by the results of gene 
sequencing studies [1,17,18]. The multifactorial nature of the 
underl3nng processes severely Hmits the intuitive interpretation 
of experimental data. Mathematical modelling is a powerful 
technique to close this gap and to provide quantitative insights 
into cell kinetics, fate determination and development of ceU 
populations. It allows a systematic study of processes not yet 
accessible by experimental procedures. Mathematical models 
have been widely applied to analyse the regulatory mechanisms 
controlling the haematopoietic system and its diseases: for 
review see [19-22] and references therein. 

The aim of this work is to investigate the impact of cell 
growth properties on the clonal selection process in acute 
leukaemias before and after treatment. We introduce math- 
ematical models of the dynamics of leukaemia, which are 
extended versions of the models proposed earlier by our 
group [23-25]. The novel ingredients of the models in this 
work are: (i) heterogeneity and multi-clonal structure of 
LSCs, (ii) different plausible feedback mechanisms and 
(iii) effects of chemotherapy. 

As the mechanism of interaction between healthy and leu- 
kaemic cell lines is not well identified, we propose two models 
(figure 1). In the first one, we assume that leukaemic cells 
depend on haematopoietic growth factors and interact with 
haematopoietic cells via competition for these factors. The 



second model is based on the assumption that autonomous 
leukaemic clones compete with haematopoietic cells for 
niches in bone marrow, which leads to increased cell death 
because of overcrowding. The latter is supported by exper- 
imental findings showing signal-independent activation of 
important cell functions [26-28] and by an increased cell 
degradation observed in leukaemic patients [29-31]. Such 
interactions have not been considered in previous models. 

The models proposed in this paper do not account for new 
mutations. Motivated by the experimental findings described 
above [5,6,8,9], we rather aim to understand which aspects of 
the dynamics of leukaemias can be explained by a selection 
process alone. It is interesting that expansion of a clone at 
relapse that could not be detected at diagnosis owing to limited 
sensitivity of methods can be misinterpreted as the occurrence 
of mutations [5]. This scenario seems to be relevant in the case 
of acute leukaemias with a short-duration treatment adminis- 
tration. Many acute leukaemias are genetically relatively 
stable in comparison to other cancers [32,33]. For this reason, 
on average many replications are necessary to acquire a new 
mutation. Consequently, it is less probable that cells acquire 
mutations during short treatment and, therefore, intrinsic 
resistance to therapy may be important, as suggested by avail- 
able evidence [5]. The latter does not hold true for long-term 
drug administration, such as imatinib treatment in the case of 
chronic leukaemias. For this reason, our work focuses on the 
acute leukaemias. For completeness of this work and to 
check how mutations might influence the model dynamics as 
it concerns results presented in this paper, we have developed 
a version of the model with mutations. The simulations of the 
model confirm our conclusions for the model with mutations. 
The model and simulations are presented in appendix A. 

Using mathematical models, we aim to identify which cell 
properties are compatible with intrinsic resistance to therapy 
and efficient expansion after treatment, and to compare 
them with the cell properties selected for before treatment. 
We perform computer simulations describing evolution of a 
multi-clonal population of leukaemic cells during disease 
development and the contribution of different clones to the 



entire cancer cell population at different time points. The het- 
erogeneity of the system is given by a certain number of 
leukaemic clones already present at the beginning of our 
observation. The models provide information on the influence 
of cell properties on the growth dynamics of the different 
clones in the presence or in the absence of chemotherapy. 
This allows us to understand how the cell properties selected 
for before treatment differ from those selected for during and 
after treatment and how treatment could be optimized to 
reduce relapses. Finally, we compare qualitatively model simu- 
lations to patients' data from clinical routine to show that the 
proposed models are consistent with clinical observations con- 
cerning the response to therapy and the time intervals between 
relapses of the disease. Details of mathematical formulation 
and parametrization of the proposed models are presented in 
appendix A. 

2. Material and methods 

2.1. Mathematical models 
2.1.1. Model assumptions 

The models used in this study are based on the models of healthy 
haematopoiesis proposed and analysed in [23,34-36] and extended 
to account for evolution of a single leukaemic clone in [25]. 

Based on the classical understanding of haematopoiesis [37], 
we assume that the system consists of an ordered sequence of 
different maturation states, so-called compartments. To describe 
time evolution of cell populations, we apply ordinary differential 
equations. The enormous amount of cells forming the haemato- 
poietic system justifies this approach [37,38]. Evolution of small 
cell population in the post-therapy period is modelled by cutting 
off the initial data which are below a minimal threshold, as it 
was, for example, proposed in [39]. 

We model time dynamics of one healthy cell lineage and an 
arbitrary number of leukaemic clones. In the description of cell 
differentiation within each cell line, we choose a two-compartment 
version of the multi-compartment system established in [23]. 
The model focuses on the maintenance of primitive cells and differ- 
entiation from undifferentiated, proliferating cells to differentiated, 
post-mitotic cells. In the case of healthy haematopoiesis, the prolifer- 
ating cells are haematopoietic stem ceUs (HSCs), haematopoietic 
progenitor cells (HPCs) and precursor cells, the post-mitotic cells 
are mature ceUs, e.g. white cells. The two-compartment architecture 
is based on a simplified description of the multi-stages differen- 
tiation process. Nevertheless, as shown in [24,35,36] models 
consisting of two compartments capture the desired d3rnamics of 
the multi-compartmental ceU population. This reduces the complex- 
ity of the differentiation process to focus on mechanisms and effects 
of competition between different cell lines. 

Each proliferating cell type is characterized by the following 
cell properties: 

— proliferation rate, describing how often a cell divides per unit 
of time; 

— fraction of self-renewal, describing the fraction of daughter 
cells returning to the compartment occupied by the mother 
cells that gave rise to them. Based on our earlier work and 
on compatibility with clinical data [23], we assume that the 
fraction of self-renewal of haematopoietic cells is regulated 
by feedback signalling; 

— death rate, describing what fraction of cells dies per unit of 
time. For simplicity, we assume that under healthy conditions 
proliferating cells do not die and post-mitotic mature blood 
cells die at a constant rate. We assume the same for leukaemic 
cells in Model 1, while in Model 2 (see below), we consider. 



additionally, cell density-dependent death rates for all bone 
marrow cell types if the marrow space is overcrowded. The 
considered marrow cell types include immature haemato- 
poietic cells and mitotic and post-mitotic leukaemic cells. 
Overcrowding is defined when marrow cell counts exceed 
the steady-state marrow cell count by two to three times. 
In this case, the death rate of post-mitotic leukaemic cells con- 
sists of their intrinsic death rate and the death triggered by 
spatial competition. 

Production of healthy blood cells is regulated by a negative 
feedback [40-42], mediated by cytokines, such as G-CSF or EPO 
[37,42,43]. If there is a shortage of blood cells of a certain type, 
the concentration of signalling molecules increases and stimulates 
expansion of precursor cells. This effect is modelled using a nega- 
tive feedback loop as proposed in [23]. Analysis and simulation 
of the model of healthy haematopoiesis, validated based on the 
clinical observations after stem cell transplantations [23,44,45], 
indicate that the regulation of the self-renewal is a more efficient 
mechanism than the regulation of the proliferation rates. Similar 
conclusions have been drawn using the models of multi-stage 
cell lineages applied to regeneration and maintenance of the 
mouse olfactory epithelium [46,47]. Therefore, in the remainder 
of this paper we assume that the regulatory mechanism is based 
on the feedback inhibition of self-renewal depending on the level 
of mature cells. 

2.1.2. Model of the healthy cell line 

We denote by p'^ the proliferation rate of mitotic haematopoietic 
cells and by the corresponding fraction of self-renewal. The 
death rate of mature blood cells is denoted by We denote 
the concentration of healthy cell types at time t by Ci(t), C2(t), cor- 
responding to mitotic and mature cells, respectively. The flux to 
mitosis at time t equals p'^{t)ci{t). During mitosis, a mother cell is 
replaced by two daughter cells. The outflux from mitosis at time t 
equals 2p''(^)ci(0, of which the fraction 2a'^{t)p''{t)ci{t) stays in 
compartment 1 (process referred to as self-renewal). The fraction 
2(1 - a{(t))p^(t)ci(t) moves to compartment 2 (process referred to 
as differentiation). 

We denote the value of the feedback signal at time t by s(f ), which 
takes values between zero and one. Self-renewal of a certain cell type 
at time t is assumed to be given as a maximal possible self-renewal of 
this ceU type multiplied by s(t). Following [23,25], we chose 
s(t) = 1/(1 + ¥c2{t)), which can be derived from cytokine kinetics 
[23]. The constant depends on the rate of extra-haematopoietic 
cytokine degradation by Hver or kidney and on the rate of cytokine 
degradation by haematopoietic cells. The latter depends on the 
densities of cytokine receptors on haematopoietic cells [45]. 

We obtain the following system of ordinary differential 
equations, where a^^^ corresponds to the maximal possible 
self-renewal of HSCs. 

^ci(0 = (2<,,s(0-l)p^ci(0, (2.1) 

^C2(0 = 2(1 - Cs(0)p^ci(f) - d^Mt) (2.2) 

and s{t) = - ^-—T. (2.3) 

The two different models proposed in this paper differ with 
respect to the interaction of leukaemic and haematopoietic ceUs. 
We consider two cases. In Model 1, leukaemic cells depend fully 
on haematopoietic cytokines, whereas in Model 2 they are totally 
independent of environmental signalling. In this sense. Models 1 
and 2 can be understood as the two opposite extremes of a conti- 
nuum. In reality, both mechanisms, competition for environmental 
signals and direct inhibition or death of haematopoietic cells, may 
contribute to impaired haematopoietic function [48]. A schematic 
of the model is given in figure 1. 



2.1.3. Model 1 

We assume that leukaemic cells depend on the same feedback 
signal as their healthy counterparts and that the post-mitotic 
leukaemic cells (blasts) decrease the supply of the factor. It 
describes a competition between healthy and leukaemic cells for 
survival signals, which results in downregulation of self-renewal. 
A schematic of the model is given in figure 1. 

To write the corresponding equations, we denote the number 
of leukaemic clones by n. As for the haematopoietic cells we 
consider mitotic and post-mitotic cell compartments for each 
leukaemic clone. Let p'' denote the proliferation rate of mitotic 
cells in leukaemic clone i and a^^^^ the corresponding maximal 
fraction of self-renewal. V>y d\> Q we denote the clearance rate 
of post-mitotic cells of clone i. Denote by l\{t) the level of mitotic 
cells of clone i and by /2(0 the level of post-mitotic cells at time t. 
These assumptions result in the following system of ordinary 
differential equations: 

d 
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(2.8) 
(2.9) 
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The expression for s{t) is a special case of s{t) = 1/ 
(1 + ¥c2{t) + Ya=i where we assume that k\ = for all L 

This simplification corresponds to the observation that the density 
of cytokine receptors is similar on cells of all leukaemic clones. For 
the major cytokine of the myeloid line, G-CSF [41], this is true for 
many patients [49]. As there is evidence that in some patients 
receptor densities may differ between different leukaemic clones 
[49], we have repeated all simulations with a randomly chosen k\ 
value for each clone, ranging from 30 to 100% of /c^. This heterogen- 
eity had no significant impact on the model results. As in many 
cases the receptor density on leukaemic cells is of the same order 
of magnitude as that on haematopoietic cells [49,50], we assume 
also k^ = ¥ for the simulation of patient examples. 



types residing in bone marrow, i.e. mitotic and post-mitotic leukae- 
mic cells as well as mitotic haematopoietic cells. For simplicity, we 
assume in Model 2 that all leukaemic cells stay in bone marrow, as 
the number of leukaemic cells exiting bone marrow is highly vari- 
able among individuals and only partially dependent on the 
leukaemia subtype [56-58] and as it is not well understood which 
mechanisms are responsible for marrow egress and high interindivi- 
dual variability. The presented results are robust with respect to this 
assumption: we repeated aU simulations for the cases that 10, 50 or 
90% of the most mature leukaemic blasts exit bone marrow. This 
has an impact on the time dynamics of marrow blast count but 
does not influence the cell properties that are selected for. 

Let d{x) be an increasing function with \m\x^ood{x) = oo. This 
function describes the death rates of bone marrow cells in depen- 
dence of bone marrow cell counts x. We assume that under 
healthy conditions there exists no cell death owing to overcrowd- 
ing. Enhanced cell death can be observed only if total bone 
marrow cellularity increases beyond the threshold level. This 
assumption is in line with bone marrow histology [59]. Therefore, 
we assume that d{x) = 0 for x < Ci, where Ci is the steady-state 
count of mitotic healthy cells. 

Assuming that the haematopoietic cell lineage is regulated as 
described above, we obtain the following system of differential 
equations: 

d 



dt 



ciit) = (2CaxS(0 - l)/ci(0 - d(x(t))c,(tl 



dt 



C2(t) = 2(l-a^^,AWci{t)-d^2C2{tl 



s{t)-- 



1 



l+¥C2(t)' 



d~i 
dt 



kit) = {Id -l)p' kit) - d{x{t))li{t\ 



(2.12) 
(2.13) 



(2.14) 



(2.15) 



it) = 2(1 - /)/ l\ it) - 4 /2 (0 - dixit))ll it), (2.16) 

; ; ; (2.17) 

plit) = i2/ -l)p^T,it)-dixit))Jlit), (2.18) 

pi it) = 2(1 - /)/ /" it) - 4 T2 (0 - dixit))~ll it) (2.19) 



and 



r(0 = ci(f) + ^?i(f) + ^?2(0- 



(2.20) 



Here ?^ denotes the mitotic cells of clone i, d their fraction of self- 
renewal and their proliferation rate. The level of post-mitotic 
cells of clone i is denoted as k- In the absence of marrow over- 
crowding, these cells die at rate d\ . 



2.1.4. Model 2 

There is evidence that in some leukaemias malignant cells show 
constitutive activation of certain signalling cascades and thus 
may become independent of external signals [26-28]. We consider 
this scenario in Model 2. In contrast to Model 1, we assume that leu- 
kaemic cells are independent of haematopoietic cytokines, whereas 
the haematopoietic cell types depend on the nonlinear feedback 
described earlier. Interaction between the healthy and cancer cell 
lines is modelled through a competition for space resulting in an 
increased cellular degradation, for example, owing to overcrowded 
bone marrow space. This is consistent with the observation of an 
increase of markers for cell death such as LDH [29-31]. Several 
mechanisms underlying this spatial competition have been pro- 
posed: (i) physical stress owing to overcrowding leads to 
extinction of cells (e.g. [51]; recently challenged by [52]), (ii) compe- 
tition for a limited niche surface expressing certain receptors 
(contact molecules) necessary for survival of the cells [53,54] and 
apoptosis if no contacts to these molecules can be established [55]. 

We model the space competition by introducing a death rate that 
increases with the number of cells in bone marrow and acts on aU cell 



2.1.5. Chemotherapy 

We focus on classical cytotoxic therapy acting on fast dividing 
cells, which is introduced to the models by adding a death rate 
proportional to the proliferation rate. The assumption is motiv- 
ated by the fact that many of the classical therapeutic agents 
used for the treatment of leukaemias act on cells in the phase 
of division or DNA replication [60]. Therefore, the rate of 
induced cell death is proportional to the number of cycling 
cells. We assume that the linear factor, denoted by /Cchemo/ is iden- 
tical for all mitotic cells. Under chemotherapy, the equation for 
mitotic haematopoietic cells in Model 1 takes the form 

^Cl(f) = (2fl'„,,s(0 - \)fc,(t) - fcchemo ' f ' C,(t). (2.21) 

Similarly, we obtain for mitotic cells of leukaemic clone i 

^^l\(t) = {TJl^^At) - 1) /m - fcchemo • p'' ■ l\(t). (1.11) 

Chemotherapy in Model 2 is introduced analogously. 
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Figure 2. Impact of growth properties on clonal selection ((a) Model 1 and (b) Model 2). The figures depict clonal selection in 50 simulated patients. Each black 
marks cellular properties of a leukaemic clone present at the beginning of the simulations in at least one patient. Each marks properties of a leukaemic clone 
contributing significantly to the leukaemic cell mass at diagnosis in at least one patient. Leukaemic cells present at diagnosis have high proliferation rates and high 
self-renewal potential. 



2.2. Simulations 

We perform numerical simulations of the models to investigate 
which leukaemic cell properties lead to survival advantage 
during evolution of leukaemogenesis and recurrence under che- 
motherapy. As explained before, the models do not account for 
additional mutations taking place during the therapy. Instead, 
we investigate evolution of a certain number of leukaemic clones 
present at a starting time point. We assume that in healthy individ- 
uals the haematopoietic cells are in a dynamic equilibrium, i.e. 
production of each cell type equals its clearance. Initial conditions 
for the computer simulations are equilibrium cell counts in the 
haematopoietic cell lineage and a small cell number for different 
leukaemic clones. We assume that the initial number of leukaemic 
clones in each patient is 50. This number is arbitrarily chosen. All 
presented simulations were repeated for different numbers of leu- 
kaemic clones (between 3 and 100), which led to comparable 
results (see the electronic supplementary material, figures SI and 
S2). We assume that primary diagnosis and diagnosis of relapse 
occur when healthy blood cell counts are decreased by 50% of 
their steady-state value. We perform simulations for 50 patients, 
i.e. 50 different sets of initial data and model parameters, with 50 
leukaemic clones per patient. The growth properties of the leukae- 
mic clones are chosen randomly within certain ranges. The choice 
of model parameters is described in appendix A. The simulations 
follow the following algorithm: 

(i) We start from healthy equilibrium in the haematopoietic line- 
age and one mitotic cell per kilogram of body weight for each 
leukaemic clone and run simulations until the number of 
healthy mature blood cells decreases by 50%. We investigate 
properties of the clones with the highest contribution to the 
total leukaemic cell mass. The clones under consideration 
are those which together constitute 80% of the total leukae- 
mic cell mass. In the following, we denote these clones as 
'significantly contributing clones'. This procedure is taken 
to reflect the sensitivity of the detection methods. In more 
than 90% of the patients, two to four clones sum up to 
more than 95% of the total leukaemic cell mass. Taking a 
threshold between 80 and 95% to define 'significantly contri- 
buting clones' has little influence on the result. Furthermore, 
more than 97% of the clones that are considered as insignifi- 
cant by this method consist of less than 1% of the leukaemic 
cell mass. This number is in agreement with the detection 
efficiency reported in the literature [61]. 

(ii) Next, we simulate chemotherapy. For simplification, we 
consider seven applications of cytotoxic drugs (one per 
day during 7 following days, corresponding to standard 
inductions). Simulations show that the number of drug 



applications has no influence on the presented qualitative 
results. As proposed in the literature [39], we assume that 
a cell population has become extinct if it consists of less 
than one cell. Initial conditions for the post-therapy period 
are obtained from cell counts after therapy where counts 
of extinct populations are set to zero. We continue simu- 
lations until mature blood cell counts decrease by 50%, 
and then assess the cell properties of the clones 
contributing to relapse. 

Calibration of the haematopoietic part of the model to clinical 
data and parameters for simulation of two patient examples can 
be found in appendix A. As in clinical routine only few key 
mutations are monitored, we choose patient examples with 
different key mutations detected at diagnosis and at relapses. 
Such data are relatively rare, therefore we focus on two patients. 
Simulations are performed using standard ODE-solvers from 
the MATLAB-software package (v. 7.8, The MathWorks, Inc., 
Natick, MA, USA) which are based on Runge-Kutta schemes. 

3. Results 

3.1. Clonality at diagnosis 

We solve the models numerically to obtain insight into the con- 
tribution of different leukaemic clones to the total leukaemic cell 
mass. Simulations indicate that at the diagnosis rarely more 
than three to four clones significantly contribute to the total leu- 
kaemic cell mass. In most cases, more than 40-50% of the total 
leukaemic cell mass originates from a single leukaemic clone. 
This finding is identical for both considered models. 

3.2. Properties of clones at diagnosis 

Simulations indicate that the clones significantly contributing 
to the leukaemic cell mass have high proliferation rates and 
high self-renewal potential (high fraction of symmetric self- 
renewing divisions). Such configuration of parameters leads 
to an efficient cell expansion. The properties of clones con- 
tributing significantly to leukaemic cell mass at diagnosis 
are depicted in figure 2. This finding is identical for both 
considered models. 

3.3. Clonality at relapse 

The clonality at relapse is comparable to the clonality at diag- 
nosis. Rarely more than three clones significantly contribute 
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Figure 3. Impact of growth properties on clonal selection. The figures depict clonal selection in 50 simulated patients. Each black marks cellular properties of a 
leukaemic clone present at the beginning of the simulations in at least one patient. Each marks properties of a leukaemic clone contributing significantly to the 
leukaemic cell mass at diagnosis in at least one patient. Grey squares mark properties of cell clones contributing significantly to relapse after chemotherapy in at 
least one patient. In comparison with leukaemic cells present at diagnosis, clones at relapse have lower proliferation rates, [a] Model 1, strong chemotherapy, 
[b] Model 1, weak chemotherapy, (c) Model 2, strong chemotherapy and (cf) Model 2, weak chemotherapy, (e) Example of the dynamics of haematopoietic 
(i) and leukaemic (ii) cells in one simulated patient. Vertical dotted lines mark primary diagnosis and relapse. Therapy is indicated by a grey rectangle. In the 
given example, primary manifestation and relapse of the disease are diagnosed when mature blood cells decreased by 50%. 



to the total leukaemic cell mass. This finding is the same for 
both considered models. 

3.4. Properties of clones at relapse 

The properties of the leukaemic clones responsible for relapse 
depend on the efficiency of chemotherapy. We run computer 
simulations for varied efficiency of chemotherapy, namely 
different death rates imposed on mitotic cell compartments. 
In the case of inefficient chemotherapy, i.e. killing rates of mito- 
tic cells being relatively small, the clones present at diagnosis are 
also responsible for relapse. These clones have high prolifer- 
ation rates and high self-renewal potential. In the case of more 
efficient chemotherapy, i.e. killing rates of mitotic cells being 



higher, the clones responsible for primary presentation differ 
from the clones responsible for relapse. Compared to the 
clones leading to primary presentation, the clones responsible 
for relapse have low proliferation rates but high self-renewal 
potential. The properties of clones contributing significantly 
to leukaemic cell mass at diagnosis and at relapse are depicted 
in figure 3. Both models lead to similar results. 

The result that slow cycling is an important selective mech- 
anism and is compatible with the finding that cells in minimal 
residual disease samples are highly quiescent [6]. It is further 
supported by the fact that addition of anthracyclines, which 
act independent of cell cycle [62], leads to improved outcome 
of relapse therapies in ALL [63]. 
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Figure 4. Time dynamics and clonal composition of subsequent relapses. The figure depicts an example of multiple relapses after chemotherapy. Relapses are 
treated using the same strategy as primary presentation, [a] Leukaemic cell counts, each shade indicates a different clone. Time between relapses 2, 3 and 4 
is shorter than remission after first treatment. This demonstrates that the selected clones are not fully responsive to the applied therapy, [b] Clonal composition 
of leukaemic cell mass at the primary diagnosis and at relapses. Charts depict the contribution of major clones to the total leukaemic cell mass. Clones responsible 
for relapse are present at very small fractions at primary diagnosis (<^5%). Relapses are triggered by the same clones but their relative contribution to the 
leukaemic cell mass change in favour of the slowly proliferating highly self-renewing cells. 



3.5. Treatment of relapse 

If the same treatment strategy as in the case of primary treat- 
ment is applied to a relapsed patient, remission time is 
significantly shorter (figure 4). Second relapse is mostly trig- 
gered by the same clones as primary relapse. With repeated 
chemotherapy, clonal composition changes in favour of the 
clones with minimal proliferation (Clone 5 in figure 4). This 
finding is in agreement with data from clinical practice in 
ALL suggesting that the clones selected for at relapse possess 
inherently reduced sensitivity to treatment [5] and may be 
also responsible for second relapse [5]. The dynamics of leu- 
kaemic cells in our model are in good agreement with data 
from clinical practice: chemotherapy is able to reduce leukae- 
mic cell load after relapses [4], nevertheless, this reduction 
does not lead to durable remission [63]. This reflects the 
worse prognosis of relapsed patients [13,63,64]. The increas- 
ing fraction of cells with reduced drug sensitivity predicted 
by the simulations explains the experimental finding that cells 
present at relapse are more resistant to chemotherapy than 
cells present at initial diagnosis [13,64]. It also shows that rep- 
etition of the same induction therapy leads to worse results in 
relapse compared with primary manifestation [63]. The selection 
of slowly cycling cells predicted by our model seems to be an 
important mechanism in AML. It was demonstrated that induc- 
tion of cell cycling enhances chemo-sensitivity of leukaemic cells 
[65] and improves patient outcome after therapy [66]. Our model 
suggests that repeated chemotherapy can lead to the selection of 
clones that are not competitive in the natural environment, i.e. 
which can be outcompeted by clones sensitive to chemotherapy 
after cessation of the treatment. 

3.6. Short-term expansion efficiency does not correlate 
with long-term self-maintenance 

If leukaemic cell behaviour depends on haematopoietic cyto- 
kines (Model 1), the current signalling environment 
influences the expansion of leukaemic clones. In this scenario. 
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Figure 5. First phase of leukaemic clone evolution: at the beginning fast 
proliferating clones with low self-renewal can dominate. They are later out- 
competed by clones with high self-renewal, which is an advantage under 
high competition for niche spaces, needed for self-renewal. If there exist 
clones with high self-renewal and high proliferation, they will dominate 
during this first phase of leukaemic evolution. Each line type corresponds 
to one leukaemic clone. Blasts are immature cells used for diagnosis of 
leukaemias. In the course of the disease, blasts accumulate and outcompete 
haematopoiesis. Blast counts greater than 5% are considered as pathological 
[67]. The simulations are based on Model 1. 

it is possible that fast proliferating cells with low self-renewal 
potential dominate the leukaemic cell mass during an initial 
phase. If, with increasing leukaemic cell mass, self-renewal 
becomes downregulated, e.g. owing to occupation of bone 
marrow niche, eventually the cell clone with the highest affi- 
nity to self -renewal survives, although its proliferation might 
be slow. An example of time evolution during an early phase 
is depicted in figure 5. 

3.7. Late relapses can originate from clones that were 
already present at diagnosis 

Simulations of Model 1 indicate that late relapses, e.g. relapses 
after more than 3 years, can originate from clones that were 



already present at diagnosis but did not significantly contrib- 
ute to the leukaemic cell mass at that time. These relapses are 
triggered by very slow proliferating cells which survive che- 
motherapy and then slowly grow. At primary diagnosis, fast 
proliferating clones dominate. The slowly proliferating clones 
are then selected by chemotherapy. This finding is able to 
explain relapses without additional mutations occurring after 
primary diagnosis. Thus, temporary risk factor exposure (e.g. 
chemicals or radiation) can also be responsible for very late 
relapses and presentations. 

3.8. Comparison of simulations to patient data 

To check whether the proposed modelling framework is 
consistent with the observed dynamics of leukaemia, we cali- 
brate the model to data of two patients with multiple relapses. 
The selected two patients showed different AML-typical 
mutations. Properties of leukaemic cells and their impairment 
owing to chemotherapy cannot be measured directly and 
the effects of specific mutations on cell d3mamics are not well 
understood. The available data include time periods between 
induction/ consolidation chemotherapy and relapse as well as 
the percentages of leukaemic blasts in the bone marrow at diag- 
nosis, follow-ups and relapse. In addition, emergence and 
subsequent elimination of leukaemia driving mutations (FLT3, 
MLL-PTD) in the bone marrow cells were precisely monitored 
using molecular biology methods [68-70]. We verify whether, 
and under which assumptions concerning the cell behaviour, 
the proposed model is compatible with clinical observations. 
This can serve as a qualitative 'proof of principle' and leads to 
hypotheses concerning changes in cell properties induced by 
the respective mutations. We assume that each mutation is 
associated with one leukaemic cell clone. We interpret differ- 
ences at diagnosis and at relapse as the result of a clonal 
selection process owing to chemotherapy and cell properties. 
For this study, we apply Model 2, as simulations over a large 
range of parameters showed that remissions shorter than 150 
days are only compatible with Model 2. 

Simulations of the evolution of leukaemic clones in the 
two patients are depicted in figures 6 and 7. The results 
show that bone marrow blast fraction can be well described 
by the model. In Patient 1, FLT3-ITD mutation of a length 
of 39 bp is detected at diagnosis. This mutation becomes 
extinct and the relapse is triggered by two different FLT3- 
ITD mutations (42 and 63 bp). This behaviour is reproduced 
in the model simulation. At diagnosis, leukaemic cell mass is 
mainly derived from one clone while at relapse two different 
clones contribute to leukaemic cell mass. 

In Patient 2, both FLT3-ITD mutation and MLL-PTD 
mutation were detected at diagnosis. The MLL-PTD mutation 
practically did not contribute to relapse. The model reflects 
this scenario. At diagnosis, two different clones contribute 
to leukaemic cell mass, one of which becomes extinct and is 
not detected at the relapse. In this patient, the clone respon- 
sible for relapse behaves similar to the HSC lineage. Thus, 
classical cytotoxic treatment would not lead to its eradication. 
This is an indication for application of new anti-leukaemic 
drugs, if feasible, or for bone marrow transplantation. 

4, Discussion 

We have examined the impact of cell properties on clonal evol- 
ution in acute leukaemias during the course of disease. 



We have considered two different mathematical models, repre- 
senting different modes of interactions between normal 
haematopoietic and leukaemic cells. In Model 1, leukaemic 
cells depend on haematopoietic cytokines, niches or other 
environmental factors. In Model 2, the leukaemic cells are inde- 
pendent of these aforementioned determinants and the only 
interaction between benign and malignant cells is owing to a 
competition for bone marrow space. 

Model simulations suggest that clones with a high pro- 
liferation rate and a high self-renewal are favoured at 
primary diagnosis. The results indicate that the number 
of clones significantly contributing to the leukaemic cell 
mass is relatively small, even if a large number of clones 
with different leukaemia driving mutations might coexist in 
the bone marrow. For example, in our simulations it was 
reduced from 50 to 2-5. This result is in agreement with 
data from recent gene sequencing studies and explains 
these data. In these studies [1,15], at most four contributing 
clones were detected in the case of AML and at most 10 in 
the case of ALL. In many patients, this number was 
even smaller. Our study implies that clonal selection owing 
to different growth characteristics is an efficient mecha- 
nism to reduce the number of clones contributing to 
leukaemic cell burden. Clones not contributing to primary 
disease manifestations might rest in a slowly proliferating 
or quiescent state and expand at relapse. Chemotherapy 
exerts a strong selective pressure on leukaemic clones, and 
thus has a considerable impact on the clonal composition 
during relapse. 

In the case of insufficient chemotherapy, the relapse can 
be triggered by the same clones as the primary disease. In 
the case of more intensive therapy regimens, relapses are 
mostly triggered by different clones than primary disease. 
This has also been concluded from experimental studies [1]. 
Our models suggest that chemotherapy selects for slowly 
proliferating clones with high self-renewal properties. 
Depending on efficiency of the therapy, it is also possible 
that clones with high proliferation and high self-renewal 
potential are responsible for relapse. 

In this study, we have focused on classical cytotoxic 
chemotherapy, mostly acting on mitotic cells. This explains 
the selection of slowly proliferating clones, among which 
those with high self-renewal potential have a competitive 
advantage, as shown in earlier studies [23-25]. High pro- 
liferation rates constitute a disadvantageous factor under 
cytotoxic treatment, as fast proliferating cells are responsive 
to even moderately intensive therapy regimens. Relapses 
due to such clones are only possible if LSCs at the same 
time have a high self-renewal potential, which is an advan- 
tageous factor for expansion and survival. Otherwise, they 
would be outcompeted by slowly proliferating cells with 
high self-renewal. Fast proliferating cells with low self- 
renewal have never been observed at relapse in our simu- 
lations. Their emergence at relapse could only be explained 
by additional mutations acquired after initial treatment. The 
selection of slowly proliferating cells may explain emergence 
of resistance in relapses. In such case, applying an identical 
therapeutic regimen to primary presentation and relapse 
has limited effects in the absence of new mutations. 

The principle of clonal competition in leukaemia evol- 
ution and the fact that resistant subclones might be 
responsible for relapse have been discussed for a long time 
[5]. Using mathematical modelling, we have provided for 
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Figure 6. Fitting of model to patient data. Different leukaemic mutations are used to distinguish between different clones, (a) The table indicates the presence and 
the absence of different leukaemic clones at different timepoints of the disease. Arrows indicate whether the respective clones increased or decreased during the 
time interval between the measurements. The depicted data are based on PCR analysis of bone marrow cells, (b) Comparison of simulated blast counts to data. Data 
are indicated as squares, (c) Evolution of leukaemic populations. Each clone is indicated by a different line type, (d) Simulated counts of healthy leucocytes. 
Chemotherapy cycles are indicated by grey rectangles. 



the first time evidence that self-renewal potential is a major 
force behind this mechanism and that cells responsible for a 
relapse show high self -renewal in nearly all cases. This find- 
ing is new and cannot be concluded from biological data 
so far. 

In appendix A, we study a model that includes occurrence 
of new mutations in addition to the selection process. In this 
scenario, the number of clones detectable at diagnosis and at 
relapse and their respective properties are practically identical 
to the scenario without mutations. This finding underlines that 
clonal selection has an important impact on the evolution of 
leukaemic cell properties. 

The exact nature of interaction between leukaemic and hae- 
matopoietic cells is not well understood. Moreover, it is well 
known that leukaemias show high interindividual heterogen- 
eity concerning symptoms and survival [67]. Therefore, it is 
possible that different mechanisms may be relevant in different 
cases. Simulation results suggest that the evolving cell proper- 
ties are robust with respect to the assumptions on the exact 



mode of interaction between haematopoietic and leukaemic 
cells and are similar in different scenarios and different 
patients. Common features of both models are: (i) relapses 
can be explained by cells that were already present at diagno- 
sis, (ii) Before therapy clonal evolution selects for cells with 
high proliferation rate and high self-renewal, (iii) Cytotoxic 
treatment selects for cells with slow proliferation and high 
self -renewal. Thus, it is possible to draw conclusions on leukae- 
mic cell properties, even if their interaction with healthy 
haematopoiesis is not known in detail. 

Nevertheless, the two proposed models exhibit some 
different dynamical properties, namely: (i) complete remis- 
sions lasting less than 150 days are only possible in Model 
2. (ii) In Model 2, it is possible that leukaemic and non-leu- 
kaemic cells coexist at ratios compatible with sufficient 
haematopoiesis for long periods, (iii) In Model 1, clones can 
temporarily expand and then be outcompeted. In Model 2, 
the clone with the fastest expansion is dominant at all times 
until treatment, (iv) In Model 2, the leukaemic cell load can 
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be reduced to a new steady state under chronic application of 
cytostatic drugs. In Model 1, expansion of leukaemic cells can 
be reduced in speed but eventually healthy haematopoiesis 
will be outcompeted. This may have application in the treat- 
ment of fast relapsing patients, as fast relapse can only be 
explained by Model 2. 

Up to now, it cannot be decided which model is more rea- 
listic. For each of the models, there exist supportive findings. 
Model 1 is supported by observations on expression of 
growth factor receptors by leukaemic cells similar to those by 
haematopoietic cells [49,50], expansion of leukaemic cells in 
the presence of cytokines in some patients [71] and dependence 
of leukaemic cell self -renewal and proliferation on chemokines 
needed for haematopoietic cell maintenance [72]. The facts sup- 
porting Model 2 are enhanced cell death in marrow samples 
[73] and increased markers for cell death/ cell lysis in serum 
[29,74], independence of leukaemic cells from important 
environmental signalling cues in the presence of some 
mutations [26] and necessity of physical contacts to marrow 
stromal cells needed for cell survival [53-55]). 



Our models support the hypothesis that processes of 
clonal selection are important mechanisms of leukaemia 
relapse, which can be responsible for expansion of different 
cell clones without a need for new mutations. A testable pre- 
diction of our models is that more sensitive methods should 
reveal larger numbers of different clones that exist but do not 
significantly contribute to the leukaemic cell mass. Another 
prediction is that cells present at relapse show mutations 
responsible for high self-renewal. 

Calibration of the models to patient data shows that 
the proposed framework is compatible with the observed 
clinical course in the considered two datasets. The predicted 
selection of slowly proliferating cells with high self-renewal 
ability is consistent with clinical observations. Our results 
may have relevance for personalized medicine. Deep sequen- 
cing techniques might provide information on the genetic 
interdependence of the clones present at diagnosis and relapse 
[1]. Our model suggests that insufficient therapy may lead 
to the presence of the same clones at diagnosis and relapse. If 
the clones present at diagnosis and relapse are not identical 



but related, i.e. they share common somatic mutations [1], 
relapse may be due to a selection process. In this case, it is prob- 
able that the clones present at relapse show a slow proliferation 
and a high self-renewal. One possible implication might be the 
application of cell-cycle independent drugs, such as those used 
in targeted therapies. 
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Appendix A 

A.I. Calibration of the haematopoietic cell lineage 

In the absence of leukaemic clones, both considered models 
reduce to the same model of the haematopoietic system. In 
steady state, this model has the following form: 

0 = (2a's-l)p' cy (Al) 
0 = 2(l-a's)p'ci-d'2C2 (A 2) 

1 



and 



Assume we know Ci and C2 . It holds 



and 
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Cl 

ci = '- 



2a' 



(A3) 

(A4) 
(A 5) 



Knowing a^, we can calculate ¥ = {2a' — l)/c2, such that 
the steady-state population size C2 is satisfied. We calibrate the 
model to the data on production of neutrophil granulocytes, 
which constitute the majority of mature white blood cells 
(50-70%). L3nTLphocytopoiesis is a complicated process involving 
l3nTiphatic organs and not only the bone marrow [75]. Therefore, 
we restrict ourselves to the myeloid line. It holds for the steady- 
state count of neutrophils [76], C2 ^ (3 - 5.8) x lOV/. We 
interpret Ci as the total amount of mitotic neutrophil precursors 
in bone marrow. Based on the data from [59], we assume that 
about 20% of bone marrow cells are mitotic precursors of neutro- 
phils (the interindividual variations are considerable). The total 
bone marrow cellularity is about 10^° cells per kilogram of 
body weight [77]. Therefore, we take 

Cl 2 X 10^ kg-^ (A 6) 

Assuming an average blood volume of 6 1 and an average 
body weight of 70 kg, we calculate a mature neutrophil count of 



C2 X 10^ kg-^ 



(A 7) 



Neutrophils have a half -life in the bloodstream, T1/2, of 
about 7h [78]. From this, we calculate 



We obtain 
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i.e. about once per 1.5 days. We know from bone marrow 
transplantation [79] that patients need about 15 days to recon- 
stitute to 5 X 10^ neutrophils per litre of blood (4 x 10^ kg~^) 
after infusion of 5 x 10^ immature cells per kilogram of body 
weight. For 

a' ^0.87 (A 10) 

this constraint is met. 



A.2. Model parameters 

For the haematopoietic branch, we chose parameters 
obtained from the calibration in §A.l. For simplicity, we 
assume = for the feedback mechanism in Model 1. We 
set the clearance rate of blasts (in the absence of effects of 
overcrowding) to d2 =0.1. This is based on the apoptotic 
indices (fraction of dying cells) reported in the literature 
which are ^0.19 ± 0.16 (19 ± 16%) [80,81]. Choosing blast 
clearance between 0.1 and 0.5 changes the speed of leukaemic 
cell accumulation but, as revealed by additional simulations, 
not the cell properties that are selected. 

We chose d{x) = iiconst. • max {0, x — Xmaxl- In histological 
images of healthy adult bone marrow, a large part of the 
bone marrow cavity consists of fat and connective tissue and 
is free of haematopoietic cells. To reflect this fact, we set 
^max ~ 2 Cl, where Ci is the steady-state count of mitotic healthy 
cells. In the simulations, dconst was set to 10~^°. This choice 
implies that if bone marrow cell counts are three times higher 
than in the steady state, the additional death rate owing to over- 
crowding is of the order of magnitude of mature cell clearance. 
The results remain unchanged qualitatively, even if values of 
^^const. vary within different orders of magnitude. 

We apply chemotherapy on 7 following days for 2 h. Differ- 
ent treatment intervals lead to comparable results. In the 
depicted simulations /Cchemo has been set to values between 
20 and 30 for less efficient therapy and to values between 40 and 
60 for efficient chemotherapy. For different choices similar 
results are obtained. The higher /Cchemo/ the stronger the selec- 
tion for high self-renewal and slow proliferation and 
the lower the probability of relapse. The lower /Cchemo/ the 
higher the probability that clones contributing to primary 
presentation are among clones contributing to relapse. For 
^chemo between 40 and 45, the obtained results are similar to 
the results obtained in experiments [1]. The parameters of the 
leukaemic clones are chosen randomly from uniform distri- 
butions, assuming that cells divide at most twice per day and 
that self-renewal is between zero and one. 

A.3. Calibration to patient examples 

Both patients were treated within clinical trials at the Univer- 
sity Hospital of Heidelberg after obtaining their written 
consent. Details on the patients' characteristics and therapeutic 
regimens can be found in the electronic supplementary 
material, table SI. Model parameters can be found in the elec- 
tronic supplementary material, tables S2 and S3. For both 
patients, the presence of specific key mutations was assessed 
in clinical routine. We chose cases where mutations get lost 
because of treatment and new mutation are detected at relapse. 
We interpret this as the result of clonal evolution. As the two 
patients harbour different mutations, their leukaemic cells 
can have different properties. Chemotherapy is modelled by 
increasing death rates for mitotic cell types during the duration 



of each cycle. For simplicity, we did not model kinetics of single 
chemotherapeutic agents. Instead, the therapy-induced death 
rates are assumed to remain constant from the first to the last 
day of each treatment cycle. In pharmacology, the exposition 
to a drug is measured using the 'area under the curve' 
(AUC). This is the integral of concentration (or drug effect) 
over time [82]. The AUC in our case is /Cchemo ' where is 
the period of drug action. The AUC over 1 day of therapy 
is similar for the single patient examples and the simulations 
in figure 3. Only myeloablative treatment before transplan- 
tation has a higher AUC. The presented results are based on 
Model 2. Model 1 is not compatible with remissions lasting 
less than 150 days. For simplicity, we count all leukaemic cell 
types as blasts. 



A.4. Model with mutations 

In the following, we describe an extension of Model 1 which 
includes mutations. There is an evidence that a preleukaemic 
HSC compartment serves as a reservoir of accumulated 
mutations [33]. This hypothesis is supported by the finding 
that some of the mutations recurrently observed in leukaemia 
already exist in the HSC compartment of a majority of leukae- 
mia patients [33]. These preleukaemic HSCs seem to behave 
similar to normal HSC [33]. The hypothesis is that a rela- 
tively small number of additional hits may transform these 
preleukaemic HSCs into LSCs. Nevertheless, details of the 
underlying dynamics are not well understood [33]. 
We make the following assumptions: 

— LSC with new properties can be generated either from 
preleukaemic HSC or from LSC owing to acquisition of 
mutations. This is in line with the current knowledge 
[7,32]. 

— For simplicity, we assume that the influx of new LSC 
from the preleukaemic compartment is constant in time. 
This assumption is made due to simplicity, because at 
the moment the dynamics of the preleukaemic compart- 
ment is not well understood [33]. We neglect mutations 
leading from normal HSC, i.e. non-preleukaemic HSC 
to LSC. 

— We assume that most mutations in LSC are acquired 
during replication of the genome and neglect other poss- 
ible origins. In line with [7], we neglect mutations leading 
to dedifferentiation of non-LSC leukaemic cells. 

Let l\ (t) be the level of LSC of clone i at time t. The flux to 
mitosis is then l\{t)p^\t). Out of mitosis, we obtain 
2a^\t)p^\t)l\(t), where a^\t) is the fraction of LSC self-renewal 
of clone i at time t. We assume that the fraction v of these cells 
is mutated, v takes into account replication errors in relevant 
genes and is assumed to be constant. The influx ai{t) of 
mutated LSCs due to new mutations occurring in clone i at 
time t is therefore 2a^\t)f (t)l[(t)v. 



We obtain the following set of equations describing 
dynamics of clone i: 

^^l\{t) = 2/{t)/{t)l\{t){l -v)- 

^/^(f) = 2(1 -/(0)/w/iw- 4/2(0 

and ai{t) = 2a^\t)/{t)l[{t)v. 

A similar system of equations has been obtained by 
Traulsen et al. [83]. As l\ is considered to be post-mitotic, 
we do not distinguish between cells that acquired a mutation 
during the divisions and those that did not. 

The influx a{t) of mutated cells at time t is given by 
0L{t) = 7 + Yl!i=i where 7 is the constant influx from 

the preleukaemic compartment and N{t) the number of 
leukaemic clones present at time t. 

We accept the rate a{t) as the rate of an inhomogeneous 
Poisson process. Poisson processes describe rare events 
[84,85], therefore they are a suitable tool to model mutations. 
It is known from probability theory that, if ti is a jump time 
of the inhomogeneous Poisson with rate X{t), then the next 
jump time T2 can be generated by solving the equation 

X{t)dt = — log(l — u) for T2. Here w is a uniformly distribu- 
ted random variable w G [0, 1] [86]. We further know that if u 
is uniformly distributed in w G [0, 1], then — log(l — u) is 
exponentially distributed with parameter 1 [85]. 

We simulate the system with mutations as follows: at time 
to = 0 we draw an exponentially distributed random number 
fi with parameter 1. We simulate the system until the time- 
point ti which fulfils f^^ \{t) dt = ri. At timepoint f 1, a 
mutation occurs that gives rise to a new LSC. This is mod- 
elled by adding to the system a new LSC clone, consisting 
of one LSC and no less primitive leukaemic cells. We 
assume that the mutation occurs in a random gene position 
therefore we assign random cell properties to the new 
clone, i.e. self-renewal and proliferation chosen randomly 
from uniform distributions (proliferation rate between 0.01 
and 0.9, self-renewal between 0.5 and 1). This choice is 
made for the sake of simplicity as details about the impact 
of mutations on cell behaviour and the underlying prob- 
ability distributions are not known [1]. We then draw 
another random number r2 and continue simulations until 
timepoint ^2 fulfilling f^^ \{t)dt = r2, etc. We start simulations 
from the equilibrium of the haematopoietic system and one 
LSC with random properties. 

The results obtained from these simulations are similar to 
the results from the model without mutations. At primary 
diagnosis, cells show high self-renewal and high prolifera- 
tion while at relapse cells show high self-renewal and 
reduced proliferation (electronic supplementary material, 
figure S3). The proliferation rates differ significantly between 
diagnosis and relapse {ip < 10~^ in Kruskal-Wallis test), 
while self -renewal does not differ significantly in 
Kruskal-Wallis test). 
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